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1.  Objective 


The  objective  of  this  research  is  to  develop  an  inversion  method  to  reconstruct  the  unknown 
atmospheric  release  parameters,  including  the  release  location  and  strength  and  associated 
concentration  field,  in  an  event  of  a  terrorist  attack.  The  inversion  method  is  based  on  the 
analysis  of  the  data  collected  from  wind,  chemical/biological  (CB)  sensors.  The  retrieved 
atmospheric  release  location  and  strength  by  this  method  should  be  the  optimal  estimations  of  the 
physical  parameters.  The  reconstructed  source  characteristics  and  concentration  field  are  useful 
information  for  force  protection  and  emergency  responses  such  as  delivering  medical  treatments, 
disinfecting  affected  areas,  and  analyzing  forensic  evidence.  The  retrieved  source  characteristics 
are  also  the  necessary  data  for  forecasting  CB  transport. 


2.  Approach 


This  source  reconstruction  system  integrates  the  backward  puff  trajectory  and  the  Bayesian 
inference  methods  to  retrieve  the  physical  parameters.  The  system  consists  of  the  data  from  CB 
detection  sensors  and  the  wind  field  from  a  diagnostic  wind  model,  the  three-dimensional  (3-D) 
wind  field  (3-DWF)  (i).  The  3-DWF  interpolates  a  limited  number  of  wind  observations  to  the 
computational  domain  using  the  mass  conservation  as  a  constraint.  A  backward  puff  trajectory 
using  the  mean  wind  field  from  the  3-DWF  is  computed  to  trace  back  to  an  approximated  area  of 
release  location.  The  improbable  release  area  is  quickly  eliminated  from  the  backward 
trajectory.  This  step  is  also  used  as  a  detection  step  to  ensure  that  the  release  location  is  in  the 
computational  domain.  After  the  backward  trajectory  computation,  the  CB  sensor  and  wind  field 
data  are  integrated  together  to  find  the  exact  release  location  and  strength  using  a  dynamical 
Bayesian  inference  theory.  The  Bayesian  theory  gives  a  maximum  likelihood  estimate  of  the 
release  parameters.  During  the  Bayesian  process,  a  forward  model  dispersion  model,  named  the 
Lagrangian  Gaussian  puff  model  (LGPM)  (2,  3),  is  applied  to  compute  the  concentration  at  each 
iterative  cycle.  The  values  of  concentration  from  CB  sensors  and  the  LGPM  model  are 
compared  and  used  to  compute  a  likelihood  function.  The  converged  state  from  the  iteration  is 
the  maximum  likelihood  estimate  of  the  source  parameters. 

The  backward  puff  trajectory  uses  the  mean  wind  field  computed  from  the  3-DWF  model  and 
parameterized  turbulence.  The  velocity  of  a  puff  centroid  at  (xp,  t)  in  a  forward  time  frame  can 
be  computed  as  u  =  dxp  /dt.  The  backward  trajectory  can  be  computed  using  a  time  coordinate 
tb=  To- t,  where  To  is  an  arbitrary  starting  time.  The  velocity  «/,  of  the  puff  in  the  backward  time 
frame  is  as  follows: 


ui,  =  dxp  /d(T o- 1)  =-  dxp  /dt. 


(1) 
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By  integrating  equation  1  backward  with  respect  to  time,  a  backward  puff  trajectory  line  can  be 
computed.  In  order  to  take  care  of  the  situation  of  an  unsteady  wind  field,  the  wind  field  needs 
to  be  frequently  updated  from  the  3-DWF  and  observations. 

The  puff  turbulence  spread  parameters,  <7f,  and  oz,  can  be  estimated  from  the  distance  between 
the  sensors  with  maximum  and  minimum  readings.  This  method  requires  that  multiple  sensors 
are  available  for  CB  detection.  In  general  CB  detection  practice,  the  sensors  are  arranged  in  an 
arc  shape  across  the  mean  wind  (a  near-circular  arrangement  for  arbitrary  wind  direction).  The 
distance  between  the  source  and  the  detection  sensor  can  be  estimated  from  formulas  from 
published  literatures  in  other  studies  (4,  5).  The  distance-spreading  relationship  is  related  to  the 
atmospheric  stability  conditions,  with  wider  spreading  in  unstable,  slower  wind  conditions,  and 
narrower  spreading  in  stable  and  stronger  wind  conditions. 


After  searching  for  the  source  area  using  the  puff  backward  trajectory  computation,  the  source 
inversion  system  applies  the  dynamical  Bayesian  inference  theory  ( 6 ,  7)  to  “fine-tune”  the 
release  location  and  find  the  release  characteristics.  Let  a  be  a  time  series  of  the  dispersion 
parameter  vector  (the  source  location  and  strength),  a  =( «/,  at ),  and  p  be  a  time  series  of 

detected  concentration  data,  P=(  pt).  The  posterior,  n(a),  a  conditional  probability  of 

system  parameters  with  respect  to  the  observational  data,  can  be  expressed  as  follows  using  the 
Bayesian  rule: 


7t((L)  = 


J9(p/«)/?(«) 

P(») 


OC  p($/a)p(a) , 


(2) 


where  p( p/a  )  is  the  likelihood  (the  conditional  probability  function  of  model  output  with  respect 
to  model  parameters),  p( a)  is  the  prior  distribution,  and  p(\\)  is  the  marginal  probability 
distribution  of  p.  The  computation  of  p( P)  is  prohibitively  expensive.  Instead  of  trying  to 
compute  it  directly,  an  alternative  approach  is  to  generate  a  collection  of  realizations  from  the 
posterior  distribution  and  use  these  samples  to  conduct  inference  (6,  7).  The  likelihood  function 
p( p/a  )  accounts  for  the  information  of  the  forward  dispersion  model  and  detection  sensor.  It 
needs  special  treatment  for  the  errors  due  to  the  sensor  threshold  limits  and  range  limit.  When 
the  concentration  is  below  the  sensor  detection  threshold,  the  reading  from  the  sensor  is  zero. 
When  the  concentration  is  exceeding  the  sensor  saturation  value,  the  instrument  reading  is  set  to 
the  saturation  level.  Following  the  treatments  of  other  investigators  (8-11),  the  likelihood 
function  p( p/a )  can  be  expressed  as  the  following  function: 


p(  p/a)  =  exp 


(logAT. -1ogC;/)2 
2  ncr2 


(3) 


where  My  and  Cy  are  the  observed  concentration  and  forward  model  produced  concentration, 
respectively,  at  location  i  and  time  j,  n  is  the  number  of  sensors  in  the  detection  network,  and  a  is 
the  error  parameter,  which  incorporates  both  the  errors  in  the  sensor  detection  and  in  the  forward 
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model.  Other  forms  of  likelihood  function  can  be  used  (7, 10).  The  function  is  designed  to 
capture  the  error  information  from  the  observations  and  the  forward  dispersion  model.  The  error 
estimate  in  Bayesian  inference  is  a  complex  research  topic;  further  research  is  required.  In  the 
current  experiment,  a  value  of  a  =  0.15  is  taken  for  the  combined  error  from  observations  and  the 
forward  model.  A  larger  value  of  a  will  yield  a  broader  final  posterior  distribution. 

The  posterior  distribution  of  nt(ai:t )  is  generated  from  the  Markov  stochastic  sampling 
procedure,  named  the  Metropolis-Hasting  method  (7, 12),  and  computed  from  equation  2.  The 
specific  Metropolis  sampling  procedure  for  our  system  can  be  described  as  follows:  The 
sampling  process  begins  with  a  source  parameter,  a  i=  (Xh  Rj),  where  X/  and  R  /  are  the  source 
location  vector  and  release  strength  at  step  1.  The  Xi  starts  at  the  previous  backward  trajectory 
obtained  region,  and  i?/  starts  at  a  reasonable  guess  of  release  rate.  In  the  i-th  iteration,  a 
candidate  state  a  =(A,-,  Rt)  is  sampled.  The  samples  of  unknown  source  parameters  ( Xj  and  ) 

are  from  a  large  set  of  possible  Xand  R.  These  source  parameters  provide  source  location  and 
release  strength  for  the  forward  dispersion  model  to  compute  the  concentration  field  at  the  sensor 
locations.  The  observed  values  and  forward  model  prediction  at  the  sensor  locations  are 
compared.  Based  on  the  comparison,  the  probability  of  source  parameter  Xi  and  Rt  at  i-th 
iteration  is  updated  using  the  likelihood  function  and  prior  probability  values  of  a  =(Xuh  Rj_i). 

If  the  comparison  of  the  Xi  and  R,  are  more  favorable  than  the  previous  (i-l)th  value,  the  sampled 
parameter  is  then  retained  for  the  next  iteration  step.  If  the  Xi  and  Rt  are  less  favorable 
parameters,  they  are  not  rejected  automatically  but  determined  by  a  random  process  with  a 
uniform  distribution.  This  ensures  that  the  sampling  process  is  not  going  to  trap  into  a  local 
optimal  value  (7,  8).  During  the  iteration,  a  candidate  source  parameter  (a* )  with  probability 
p( a,  a  )  is  computed  in  the  follow  equation: 

p(a, a  )  =  mm(f| \  1) ,  (4) 

«=i  A,(«) 

where  n  is  the  number  of  sensors  in  the  network  and  ^(aj  and  nn  (a)  are  computed  using 

equation  2.  Typically,  four  Markov  chains  are  used  in  the  stochastic  sampling  process.  The 
convergence  is  attained  when  the  ratio  of  variance  within  chain  to  variance  between  chains  is 
approaching  unity  (6).  The  final  posterior  distribution  produces  the  most  likely  source 
parameters  X  and  R.  The  final  forward  concentration  prediction  is  the  reconstructed 
concentration  field. 


3.  Results 


Before  applying  the  inversion  method  to  a  complex  urban  environment,  an  experiment  with  a 
simple  situation  was  performed.  In  this  example,  a  release  source  reconstruction  in  an  idealized 
suburban  area  was  considered.  The  computation  domain  consisted  of  five  buildings  over  a 
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gentle  terrain  (figure  1).  The  chemical  detection  sensors  were  arranged  at  ground  surface  in  a 
circle  enclosing  the  five  buildings.  The  test  was  based  on  a  simulated  “truth”  of  an  atmospheric 
release  at  the  ground  surface.  The  wind  field  was  initialized  with  several  logarithmic  profiles. 
Applying  the  data  of  the  terrain,  the  building  height,  and  the  initial  wind  profiles,  the  wind  field 
was  computed  using  the  3-DWF  model.  The  wind  field  remained  in  a  steady  condition  during 
the  20-min  release  period.  A  neutrally  buoyant  gas  was  released  in  an  upwind  area  of  the  five 
buildings.  The  concentration  field  was  simulated  with  a  Lagrangian  stochastic  particle  model 
(LSPM)  (13-15),  which  was  to  be  compared  with  the  reconstructed  field.  The  LSPM-type 
model  is  usually  considered  a  more  advanced  and  accurate  approach  for  the  air  pollutant 
dispersion  simulation  (13,  14).  Unlike  the  LGPM  model,  a  large  number  of  particles  in  LSPM 
are  released  and  the  turbulence  effects  are  simulated  with  a  random  walk  model.  However,  a 
much  longer  computation  time  (as  much  as  10x  more  time  compared  with  the  LGPM)  is  required 
for  a  LSPM  model  simulation.  Details  about  the  LSPM  can  be  found  in  Wilson  and  Sawford 
(13)  and  Thomson  (14). 


(a) 


Y(km)  0  o  X<km> 


Figure  1.  A  3-D  diagram  for  the  simulation  domain. 

Assume  that  in  a  real  situation,  the  release  source  location  and  release  strength  are  not  known, 
only  the  concentration  time  series  sampled  (figure  2)  from  6  out  of  the  8  sensors  (in  this 
simulated  case).  The  wind  field  condition  is  monitored  continuously.  From  the  inspection  of  the 
time  series  signal  (figure  2)  and  the  wind  field  (figure  3),  the  release  is  continuous  in  a  rather 
steady  wind  condition  since  the  concentration  signals  reached  the  constant  values  for  each  sensor 
after  about  100  s.  The  inverse  system  uses  the  information  from  the  detection  sensors  and 
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Figure  2.  The  time  series  of  concentration  detected  by 
sensors;  sensors  1  and  8  have  blank  readings. 


Figure  3.  The  wind  field  computed  from  3-DWF  at 
4  m.  The  dashed  line  denotes  the  puff 
backward  trajectory,  and  the  puff  spreading 
parameter  is  represented  by  circle  radius. 
The  red  dot  is  the  starting  point  of  a  Markov 
chain,  and  the  red  square  is  the  release 
location  by  the  inversion  method. 
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the  wind  field  to  reconstruct  the  source  location  and  the  release  strength.  The  first  step  is  to  do  a 
backward  trajectory  computation,  starting  at  the  sensor  location  with  the  maximum  reading. 

This  sensor  location  is  considered  as  an  approximated  plume  centroid  location.  The  starting 
plume  parameter,  a/„  is  approximated  from  the  line  which  is  perpendicular  to  the  wind  at  the 
sensor  location.  The  oz  parameter  is  taken  as  0.75a/,  to  start  the  backward  trajectory  computation 
for  a  neutral  atmospheric  stability  condition.  The  plume  source  is  traced  back  using  the 
backward  trajectory  shown  in  figure  3,  where  the  radiuses  of  the  circles  represent  the  Oh  values 
of  the  puff.  The  Oh  and  az  values  are  computed  using  the  formula  from  previous  studies  (5,  6). 
The  puff  centroid  backward  trajectory  is  computed  by  solving  equation  1  using  a  second-order 
Runge-Kutta  method.  Since  the  purpose  of  the  backward  trajectory  step  is  to  reduce  the 
computation  in  the  entire  source  inversion  system  by  tracing  back  to  the  approximated  source 
location,  the  trajectory  computation  is  stopped  at  the  location  with  the  <7/,~  0.1  km.  At  this  time, 
the  Markov  chain  Stochastic  sampling  process  takes  over  to  find  the  maximum  likelihood 
estimate  of  release  location  and  characteristics. 

Four  Markov  chains  started  at  the  area  where  the  puff  backward  trajectory  was  traced.  The 
Markov  chains  sampled  source  locations  and  release  strength,  following  the  rules  described  in 
the  Metropolis-Hasting  method  (7,  12).  For  clarity,  figure  3  only  shows  the  sampling  locations 
from  one  Markov  chain,  where  the  starting  point  is  marked  by  a  red  circle.  The  converging  point 
is  denoted  by  a  small  red  square.  The  sampling  starting  point  is  proven  to  be  near  the  release 
point  with  the  help  of  backward  trajectory.  If  the  reconstruction  system  was  not  preceded  with  a 
backward  trajectory  computation,  the  Markov  chain  stochastic  sampling  would  take  much 
longer.  The  sampling  starting  location  could  have  started  at  a  location  far  away  from  the  release 
location,  taking  much  more  iterations  to  converge  to  the  release  point.  The  system  has  been 
tested  without  using  the  backward  trajectory  procedure  starting  at  random  locations  in  the 
computational  domain.  Retrieving  the  release  location  would  take  as  much  as  4  hr  of  CPU  time 
(~3x  the  computation  time  with  default  inversion  system)  on  a  3-GH  workstation  to  converge  to 
the  release  point.  At  the  same  stochastic  sampling  process,  the  release  strength  was  statistically 
retrieved  in  the  Bayesian  procedure.  The  sampling  started  at  the  release  strength  of  1  g/s,  carried 
out  with  the  0. 1  increment  of  release  strength  for  a  sampling  iteration.  The  resulting  release 
strength  sampling  is  shown  in  a  probability  distribution  in  figure  4.  It  indicated  that  the  largest 
posterior  probability  is  approximately  corresponding  to  the  “true”  release  strength  in  the  base 
simulation.  The  maximum  likelihood  estimate  of  the  release  strength  is  slightly  greater  than 
“true”  release  strength  of  3  g/s.  Using  the  reconstructed  release  location  and  release  strength,  a 
concentration  field  is  computed  as  shown  in  figure  5.  The  comparison  with  the  base  simulation 
(figure  6)  indicated  that  the  release  was  satisfactory  in  this  simple  situation. 
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Figure  4.  The  sampling  frequency  for  the  release  strength. 
The  red  line  is  the  “true”  release  strength. 


Figure  5.  The  reconstructed  plume  using  the  inversion  method. 
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Figure  6.  The  simulated  “truth”  plume  using  the  LSPM  in 
3-DWF  computed  wind  field. 


4.  Conclusions 


A  source  inversion  method  for  reconstruction  of  the  surface  atmospheric  release  source  location 
and  strength  has  been  developed  in  this  research.  The  inversion  system  consists  of  a  backward 
puff  trajectory  computation  and  a  Bayesian  inference  process  via  stochastic  sampling.  The 
backward  puff  trajectory  computation  is  used  to  trace  back  the  approximated  area  of  release 
location.  The  approximated  area  from  the  backward  trajectory  serves  as  a  starting  point  for  a 
Bayesian  inference  process.  The  Bayesian  inference  process  is  applied  to  search  the  release 
strength  and  the  release  location.  A  preliminary  test  with  an  idealized  continuous,  point  source 
release  case  indicated  that  the  method  gives  satisfactory  inversion  results  in  terms  of  the  release 
location  and  strength.  This  integrated  methodology  combined  backward  trajectory  and  Bayesian 
inference  reduces  about  66%  of  the  computation  time  compared  with  the  method  using  the 
Bayesian  inference  only  in  a  test  case.  Obviously,  this  is  a  preliminary  exploration  of  the 
inversion  method  for  a  very  simple  situation.  For  real  complex  urban  conditions,  the  following 
difficult  issues  need  to  be  resolved:  (1)  the  complex  turbulence  characteristic  in  the  urban 
environment,  (2)  the  backward  trajectory  computation  in  an  urban  environment,  (3)  the 
computation  speed  of  the  inversion  system,  and  (4)  the  forward  modeling  of  the  plum  or  puff  in 
different  atmospheric  stability  conditions. 
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